knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")
library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE #
if (get_bayesfactor) {
stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}
options(
dplyr.print_max = 100,
brms.backend = 'cmdstan',
brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
brms.file_refit = 'on_change',
#brms.file_refit = 'always',
error = function() {
beepr::beep(sound = 5)
if (shutdown) {
system("shutdown /s /t 180")
quit(save = "no", status = 1)
}
}
, es.use_symbols = TRUE
)
####################### Model parameters #######################
iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000
# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'
################################################################
suffix = paste0('_AsPreregisteredAPIM_', as.character(iterations))df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df
df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]Constructing scales Re-coding pusing reshaping data (4field) centering data within and between
Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_persuasion", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
## Start sampling
if (check_models) {
check_brms(pa_sub_persuasion, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.24 [2.14, 2.33] 1.49 0.45
## is_B 1.06 [1.03, 1.10] 1.03 0.95
## is_A:persuasion_self_cw 1.15 [1.12, 1.19] 1.07 0.87
## is_B:persuasion_self_cw 1.14 [1.11, 1.18] 1.07 0.88
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 2.71 [2.60, 2.84] 1.65 0.37
## is_A:persuasion_self_cb 2.86 [2.74, 2.99] 1.69 0.35
## Tolerance 95% CI
## [0.43, 0.47]
## [0.91, 0.97]
## [0.84, 0.89]
## [0.85, 0.90]
## [0.92, 0.98]
## [0.35, 0.38]
## [0.33, 0.37]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 94%
## lognormal 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 78%
## beta-binomial 16%
## lognormal 3%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 7, observations = 3742, p-value =
## 0.84
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0005344735 0.0029396045
## sample estimates:
## outlier frequency (expected: 0.001574024585783 )
## 0.001870657
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_persuasion$data$pa_sub[pa_sub_persuasion$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_persuasion,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.74 | 0.14 | [ 0.51, 1.07] | 0.944 | [0.84, 1.20] | 0.257 | 0.260 | Moderate Evidence for Null | 1.000 | 17468 | 24772 | 46.58*** | 3.50 | [40.07, 53.85] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 13642 | 22685 |
| is_A:day | 0.97 | 0.18 | [ 0.67, 1.40] | 0.570 | [0.84, 1.20] | 0.660 | 0.095 | Strong Evidence for Null | 1.000 | 44612 | 29960 | 1.09 | 0.09 | [ 0.92, 1.29] | 0.845 | [0.92, 1.08] | 0.436 | 0.012 | Very Strong Evidence for Null | 1.000 | 46827 | 29347 |
| is_A:persuasion_partner_cb | 0.75 | 0.30 | [ 0.34, 1.64] | 0.767 | [0.84, 1.20] | 0.271 | 0.254 | Moderate Evidence for Null | 1.000 | 14088 | 21142 | 1.02 | 0.17 | [ 0.74, 1.42] | 0.561 | [0.92, 1.08] | 0.365 | 0.018 | Very Strong Evidence for Null | 1.000 | 11377 | 19396 |
| is_A:persuasion_partner_cw | 1.58*** | 0.16 | [ 1.30, 1.98] | 1.000 | [0.84, 1.20] | 0.004 | >100 | Overwhelming Evidence | 1.000 | 23821 | 26789 | 1.02 | 0.03 | [ 0.96, 1.08] | 0.712 | [0.92, 1.08] | 0.980 | 0.008 | Very Strong Evidence for Null | 1.000 | 20463 | 25811 |
| is_A:persuasion_self_cb | 0.81 | 0.36 | [ 0.33, 1.97] | 0.683 | [0.84, 1.20] | 0.281 | 0.251 | Moderate Evidence for Null | 1.000 | 13340 | 20895 | 1.21 | 0.21 | [ 0.84, 1.72] | 0.853 | [0.92, 1.08] | 0.199 | 0.029 | Very Strong Evidence for Null | 1.001 | 10250 | 19030 |
| is_A:persuasion_self_cw | 1.67*** | 0.16 | [ 1.41, 2.10] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 22586 | 23823 | 1.01 | 0.04 | [ 0.94, 1.10] | 0.643 | [0.92, 1.08] | 0.936 | 0.011 | Very Strong Evidence for Null | 1.000 | 15189 | 22445 |
| is_B | 0.92 | 0.18 | [ 0.63, 1.37] | 0.656 | [0.84, 1.20] | 0.597 | 0.084 | Strong Evidence for Null | 1.000 | 16516 | 22811 | 49.76*** | 3.97 | [42.49, 58.38] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 11733 | 19760 |
| is_B:day | 0.87 | 0.16 | [ 0.61, 1.25] | 0.777 | [0.84, 1.20] | 0.543 | 0.123 | Moderate Evidence for Null | 1.000 | 41627 | 29767 | 0.88 | 0.07 | [ 0.75, 1.03] | 0.943 | [0.92, 1.08] | 0.257 | 0.024 | Very Strong Evidence for Null | 1.000 | 50533 | 30069 |
| is_B:persuasion_partner_cb | 1.01 | 0.48 | [ 0.39, 2.62] | 0.510 | [0.84, 1.20] | 0.293 | 0.235 | Moderate Evidence for Null | 1.000 | 12851 | 19993 | 1.03 | 0.19 | [ 0.72, 1.48] | 0.568 | [0.92, 1.08] | 0.330 | 0.017 | Very Strong Evidence for Null | 1.000 | 11737 | 19703 |
| is_B:persuasion_partner_cw | 1.49*** | 0.14 | [ 1.24, 1.85] | 1.000 | [0.84, 1.20] | 0.010 | 86.845 | Very Strong Evidence | 1.000 | 26845 | 24277 | 1.03 | 0.03 | [ 0.97, 1.10] | 0.858 | [0.92, 1.08] | 0.935 | 0.014 | Very Strong Evidence for Null | 1.000 | 20439 | 27855 |
| is_B:persuasion_self_cb | 0.70 | 0.29 | [ 0.30, 1.62] | 0.805 | [0.84, 1.20] | 0.233 | 0.304 | Weak Evidence for Null | 1.000 | 13544 | 21173 | 1.13 | 0.18 | [ 0.82, 1.57] | 0.767 | [0.92, 1.08] | 0.289 | 0.023 | Very Strong Evidence for Null | 1.000 | 12920 | 20807 |
| is_B:persuasion_self_cw | 1.68*** | 0.17 | [ 1.39, 2.10] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 22677 | 24944 | 1.04 | 0.03 | [ 0.98, 1.10] | 0.889 | [0.92, 1.08] | 0.930 | 0.016 | Very Strong Evidence for Null | 1.000 | 17302 | 24439 |
| sd(is_A) | 0.94 | 0.12 | [0.73, 1.23] | NA | NA | NA | NA | NA | 1.000 | 20183 | 25071 | 0.33 | 0.05 | [0.25, 0.44] | NA | NA | NA | NA | NA | 1.000 | 13315 | 23621 |
| sd(is_A:persuasion_partner_cw) | 0.41 | 0.12 | [0.21, 0.69] | NA | NA | NA | NA | NA | 1.000 | 18656 | 21846 | 0.06 | 0.03 | [0.01, 0.13] | NA | NA | NA | NA | NA | 1.000 | 12477 | 10932 |
| sd(is_A:persuasion_self_cw) | 0.29 | 0.13 | [0.05, 0.58] | NA | NA | NA | NA | NA | 1.000 | 9601 | 8514 | 0.17 | 0.03 | [0.11, 0.24] | NA | NA | NA | NA | NA | 1.000 | 20343 | 28866 |
| sd(is_B) | 1.01 | 0.13 | [0.79, 1.32] | NA | NA | NA | NA | NA | 1.000 | 18297 | 24221 | 0.39 | 0.05 | [0.30, 0.50] | NA | NA | NA | NA | NA | 1.000 | 13771 | 21334 |
| sd(is_B:persuasion_partner_cw) | 0.33 | 0.14 | [0.06, 0.63] | NA | NA | NA | NA | NA | 1.000 | 10537 | 8117 | 0.10 | 0.03 | [0.05, 0.16] | NA | NA | NA | NA | NA | 1.000 | 20368 | 18826 |
| sd(is_B:persuasion_self_cw) | 0.38 | 0.10 | [0.19, 0.62] | NA | NA | NA | NA | NA | 1.000 | 21344 | 22705 | 0.09 | 0.03 | [0.05, 0.15] | NA | NA | NA | NA | NA | 1.000 | 18840 | 17891 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.67 | 0.01 | [0.65, 0.70] | NA | NA | NA | NA | NA | 1.000 | 47834 | 30606 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pressure", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
## Start sampling
if (check_models) {
check_brms(pa_sub_pressure, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_pressure, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.23 [2.14, 2.33] 1.49 0.45
## is_B 1.05 [1.02, 1.09] 1.02 0.96
## is_A:pressure_self_cw 1.01 [1.00, 1.14] 1.01 0.99
## is_B:pressure_self_cw 1.01 [1.00, 1.12] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 1.15] 1.01 0.99
## is_B:pressure_partner_cw 3.49 [3.33, 3.66] 1.87 0.29
## is_A:pressure_self_cb 4.30 [4.10, 4.52] 2.07 0.23
## Tolerance 95% CI
## [0.43, 0.47]
## [0.92, 0.98]
## [0.88, 1.00]
## [0.89, 1.00]
## [0.87, 1.00]
## [0.27, 0.30]
## [0.22, 0.24]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 94%
## lognormal 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 78%
## beta-binomial 16%
## lognormal 3%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 11 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 6, observations = 3736, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0005353319 0.0030848501
## sample estimates:
## outlier frequency (expected: 0.00157922912205567 )
## 0.001605996
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pressure$data$pa_sub[pa_sub_pressure$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_pressure,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.87 | 0.15 | [ 0.61, 1.23] | 0.783 | [0.84, 1.20] | 0.555 | 0.096 | Strong Evidence for Null | 1.000 | 14549 | 20069 | 48.25*** | 3.51 | [41.79, 55.67] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 10797 | 19767 |
| is_A:day | 0.67* | 0.12 | [ 0.48, 0.95] | 0.989 | [0.84, 1.20] | 0.107 | 1.172 | Weak Evidence | 1.000 | 37555 | 29688 | 1.09 | 0.09 | [ 0.92, 1.28] | 0.835 | [0.92, 1.08] | 0.452 | 0.011 | Very Strong Evidence for Null | 1.000 | 48844 | 30086 |
| is_A:pressure_partner_cb | 0.63 | 0.39 | [ 0.18, 2.10] | 0.772 | [0.84, 1.20] | 0.174 | 0.412 | Weak Evidence for Null | 1.000 | 14001 | 22562 | 0.94 | 0.26 | [ 0.54, 1.64] | 0.582 | [0.92, 1.08] | 0.218 | 0.016 | Very Strong Evidence for Null | 1.000 | 12801 | 20013 |
| is_A:pressure_partner_cw | 2.44*** | 0.57 | [ 1.57, 4.32] | 1.000 | [0.84, 1.20] | 0.001 | >100 | Overwhelming Evidence | 1.000 | 27317 | 18821 | 0.93 | 0.05 | [ 0.84, 1.05] | 0.880 | [0.92, 1.08] | 0.567 | 0.012 | Very Strong Evidence for Null | 1.000 | 29892 | 25384 |
| is_A:pressure_self_cb | 0.57 | 0.44 | [ 0.13, 2.62] | 0.766 | [0.84, 1.20] | 0.143 | 0.502 | Weak Evidence for Null | 1.000 | 14022 | 21390 | 1.56 | 0.53 | [ 0.80, 3.08] | 0.908 | [0.92, 1.08] | 0.075 | 0.038 | Strong Evidence for Null | 1.000 | 14697 | 21985 |
| is_A:pressure_self_cw | 1.86* | 0.62 | [ 1.02, 4.68] | 0.978 | [0.84, 1.20] | 0.060 | 1.518 | Weak Evidence | 1.000 | 19579 | 18590 | 0.85* | 0.06 | [ 0.72, 1.00] | 0.978 | [0.92, 1.08] | 0.136 | 0.082 | Strong Evidence for Null | 1.000 | 23058 | 21388 |
| is_B | 1.14 | 0.21 | [ 0.80, 1.63] | 0.763 | [0.84, 1.20] | 0.560 | 0.090 | Strong Evidence for Null | 1.000 | 13997 | 21736 | 52.20*** | 4.30 | [44.46, 61.41] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 9866 | 16264 |
| is_B:day | 0.56*** | 0.10 | [ 0.40, 0.79] | 1.000 | [0.84, 1.20] | 0.011 | 20.363 | Strong Evidence | 1.000 | 36258 | 29993 | 0.83* | 0.07 | [ 0.71, 0.97] | 0.991 | [0.92, 1.08] | 0.085 | 0.109 | Moderate Evidence for Null | 1.000 | 48596 | 30206 |
| is_B:pressure_partner_cb | 1.13 | 0.88 | [ 0.24, 5.28] | 0.562 | [0.84, 1.20] | 0.183 | 0.393 | Weak Evidence for Null | 1.000 | 12957 | 20347 | 1.10 | 0.43 | [ 0.49, 2.47] | 0.593 | [0.92, 1.08] | 0.155 | 0.018 | Very Strong Evidence for Null | 1.001 | 12994 | 19397 |
| is_B:pressure_partner_cw | 3.28* | 2.01 | [ 1.24, 16.35] | 0.990 | [0.84, 1.20] | 0.017 | 4.367 | Moderate Evidence | 1.000 | 17196 | 17609 | 0.99 | 0.07 | [ 0.82, 1.13] | 0.575 | [0.92, 1.08] | 0.709 | 0.008 | Very Strong Evidence for Null | 1.000 | 21016 | 19395 |
| is_B:pressure_self_cb | 0.48 | 0.30 | [ 0.14, 1.67] | 0.882 | [0.84, 1.20] | 0.111 | 0.645 | Weak Evidence for Null | 1.000 | 12815 | 19815 | 1.08 | 0.34 | [ 0.57, 2.05] | 0.593 | [0.92, 1.08] | 0.194 | 0.019 | Very Strong Evidence for Null | 1.001 | 12751 | 19790 |
| is_B:pressure_self_cw | 1.29 | 0.34 | [ 0.71, 2.66] | 0.832 | [0.84, 1.20] | 0.326 | 0.227 | Moderate Evidence for Null | 1.000 | 20759 | 17702 | 0.93 | 0.06 | [ 0.82, 1.08] | 0.842 | [0.92, 1.08] | 0.545 | 0.012 | Very Strong Evidence for Null | 1.000 | 23751 | 24390 |
| sd(is_A) | 0.88 | 0.12 | [0.68, 1.15] | NA | NA | NA | NA | NA | 1.000 | 17349 | 24619 | 0.32 | 0.05 | [0.24, 0.44] | NA | NA | NA | NA | NA | 1.000 | 12153 | 20889 |
| sd(is_A:pressure_partner_cw) | 0.31 | 0.28 | [0.01, 1.20] | NA | NA | NA | NA | NA | 1.000 | 15331 | 16422 | 0.06 | 0.05 | [0.00, 0.21] | NA | NA | NA | NA | NA | 1.000 | 17943 | 17235 |
| sd(is_A:pressure_self_cw) | 0.72 | 0.53 | [0.06, 2.40] | NA | NA | NA | NA | NA | 1.000 | 11421 | 12927 | 0.10 | 0.09 | [0.00, 0.38] | NA | NA | NA | NA | NA | 1.000 | 12722 | 15622 |
| sd(is_B) | 0.91 | 0.12 | [0.70, 1.19] | NA | NA | NA | NA | NA | 1.000 | 18280 | 24673 | 0.40 | 0.05 | [0.31, 0.53] | NA | NA | NA | NA | NA | 1.000 | 11041 | 18555 |
| sd(is_B:pressure_partner_cw) | 1.65 | 0.85 | [0.31, 3.83] | NA | NA | NA | NA | NA | 1.000 | 12805 | 11132 | 0.09 | 0.08 | [0.00, 0.37] | NA | NA | NA | NA | NA | 1.000 | 13239 | 14438 |
| sd(is_B:pressure_self_cw) | 0.74 | 0.49 | [0.06, 2.12] | NA | NA | NA | NA | NA | 1.000 | 10281 | 10816 | 0.09 | 0.07 | [0.00, 0.28] | NA | NA | NA | NA | NA | 1.000 | 14609 | 13572 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.69 | 0.01 | [0.67, 0.72] | NA | NA | NA | NA | NA | 1.000 | 43470 | 29205 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
## Start sampling
if (check_models) {
check_brms(pa_sub_pushing, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_pushing, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.75 [1.68, 1.84] 1.32 0.57
## is_B 1.08 [1.05, 1.12] 1.04 0.93
## barriers_self_cw 1.04 [1.02, 1.10] 1.02 0.96
## barriers_self_cb 1.04 [1.01, 1.09] 1.02 0.97
## is_A:pushing_self_cw 1.07 [1.04, 1.12] 1.04 0.93
## is_B:pushing_self_cw 1.22 [1.18, 1.28] 1.11 0.82
## is_A:pushing_partner_cw 1.24 [1.20, 1.30] 1.12 0.80
## is_B:pushing_partner_cw 1.19 [1.15, 1.24] 1.09 0.84
## is_A:pushing_self_cb 1.27 [1.23, 1.33] 1.13 0.79
## Tolerance 95% CI
## [0.54, 0.60]
## [0.89, 0.95]
## [0.91, 0.98]
## [0.92, 0.99]
## [0.89, 0.96]
## [0.78, 0.85]
## [0.77, 0.83]
## [0.81, 0.87]
## [0.75, 0.82]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 91%
## lognormal 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 84%
## beta-binomial 9%
## lognormal 6%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 12 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 5, observations = 1776, p-value =
## 0.36
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00423705
## sample estimates:
## outlier frequency (expected: 0.00161036036036036 )
## 0.002815315
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pushing$data$pa_sub[pa_sub_pushing$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_pushing,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| barriers_self_cb | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.89 | 0.13 | [ 0.68, 1.18] | 0.791 | [0.93, 1.08] | 0.305 | 0.026 | Very Strong Evidence for Null | 1.000 | 15389 | 22099 |
| barriers_self_cw | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 1.06 | 0.04 | [ 0.99, 1.13] | 0.943 | [0.93, 1.08] | 0.739 | 0.032 | Strong Evidence for Null | 1.000 | 68293 | 29503 |
| is_A | 2.63*** | 0.70 | [ 1.56, 4.51] | 1.000 | [0.84, 1.20] | 0.002 | 48.662 | Very Strong Evidence | 1.000 | 29400 | 28949 | 53.39*** | 4.54 | [45.24, 63.19] | 1.000 | [0.93, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 18788 | 26280 |
| is_A:day | 0.82 | 0.24 | [ 0.46, 1.45] | 0.751 | [0.84, 1.20] | 0.369 | 0.187 | Moderate Evidence for Null | 1.000 | 58113 | 28964 | 1.08 | 0.10 | [ 0.90, 1.30] | 0.793 | [0.93, 1.08] | 0.441 | 0.011 | Very Strong Evidence for Null | 1.000 | 66370 | 30302 |
| is_A:pushing_partner_cb | 0.91 | 0.88 | [ 0.13, 6.16] | 0.540 | [0.84, 1.20] | 0.147 | 0.491 | Weak Evidence for Null | 1.000 | 30887 | 30885 | 1.10 | 0.40 | [ 0.52, 2.26] | 0.597 | [0.93, 1.08] | 0.159 | 0.020 | Very Strong Evidence for Null | 1.000 | 15843 | 23109 |
| is_A:pushing_partner_cw | 1.95** | 0.45 | [ 1.30, 3.48] | 0.999 | [0.84, 1.20] | 0.009 | 14.509 | Strong Evidence | 1.000 | 32690 | 26682 | 0.96 | 0.03 | [ 0.89, 1.03] | 0.886 | [0.93, 1.08] | 0.837 | 0.017 | Very Strong Evidence for Null | 1.000 | 47103 | 27861 |
| is_A:pushing_self_cb | 0.39 | 0.35 | [ 0.07, 2.31] | 0.852 | [0.84, 1.20] | 0.093 | 0.769 | Weak Evidence for Null | 1.000 | 28601 | 30135 | 1.29 | 0.42 | [ 0.67, 2.49] | 0.785 | [0.93, 1.08] | 0.136 | 0.025 | Very Strong Evidence for Null | 1.000 | 18803 | 26620 |
| is_A:pushing_self_cw | 1.91 | 0.76 | [ 0.96, 5.60] | 0.967 | [0.84, 1.20] | 0.081 | 0.998 | Weak Evidence for Null | 1.000 | 25693 | 21815 | 1.00 | 0.06 | [ 0.89, 1.12] | 0.523 | [0.93, 1.08] | 0.814 | 0.012 | Very Strong Evidence for Null | 1.000 | 24382 | 27449 |
| is_B | 4.02*** | 1.18 | [ 2.24, 7.28] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 26916 | 28674 | 52.54*** | 4.45 | [44.44, 62.18] | 1.000 | [0.93, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 20116 | 25664 |
| is_B:day | 0.65 | 0.19 | [ 0.36, 1.15] | 0.929 | [0.84, 1.20] | 0.177 | 0.427 | Weak Evidence for Null | 1.000 | 63463 | 30099 | 0.94 | 0.09 | [ 0.79, 1.12] | 0.746 | [0.93, 1.08] | 0.507 | 0.010 | Very Strong Evidence for Null | 1.000 | 74905 | 30566 |
| is_B:pushing_partner_cb | 0.81 | 0.81 | [ 0.11, 5.67] | 0.582 | [0.84, 1.20] | 0.140 | 0.508 | Weak Evidence for Null | 1.000 | 29944 | 29697 | 1.04 | 0.34 | [ 0.53, 2.01] | 0.543 | [0.93, 1.08] | 0.187 | 0.020 | Very Strong Evidence for Null | 1.000 | 19801 | 25578 |
| is_B:pushing_partner_cw | 2.28* | 1.11 | [ 1.04, 8.55] | 0.979 | [0.84, 1.20] | 0.044 | 1.868 | Weak Evidence | 1.000 | 21839 | 21728 | 1.02 | 0.06 | [ 0.91, 1.15] | 0.658 | [0.93, 1.08] | 0.797 | 0.012 | Very Strong Evidence for Null | 1.000 | 27666 | 30385 |
| is_B:pushing_self_cb | 0.41 | 0.42 | [ 0.05, 3.18] | 0.810 | [0.84, 1.20] | 0.096 | 0.757 | Weak Evidence for Null | 1.000 | 31398 | 29450 | 1.50 | 0.55 | [ 0.71, 3.13] | 0.859 | [0.93, 1.08] | 0.092 | 0.038 | Strong Evidence for Null | 1.000 | 17644 | 25436 |
| is_B:pushing_self_cw | 1.16 | 0.25 | [ 0.73, 1.91] | 0.753 | [0.84, 1.20] | 0.491 | 0.137 | Moderate Evidence for Null | 1.000 | 29855 | 26726 | 0.95 | 0.03 | [ 0.88, 1.01] | 0.939 | [0.93, 1.08] | 0.726 | 0.030 | Very Strong Evidence for Null | 1.000 | 43301 | 31543 |
| sd(is_A) | 1.18 | 0.20 | [0.85, 1.64] | NA | NA | NA | NA | NA | 1.000 | 27290 | 31419 | 0.36 | 0.06 | [0.26, 0.50] | NA | NA | NA | NA | NA | 1.000 | 16350 | 24786 |
| sd(is_A:pushing_partner_cw) | 0.65 | 0.33 | [0.09, 1.45] | NA | NA | NA | NA | NA | 1.000 | 14902 | 13091 | 0.04 | 0.04 | [0.00, 0.14] | NA | NA | NA | NA | NA | 1.000 | 20863 | 21692 |
| sd(is_A:pushing_self_cw) | 1.51 | 0.59 | [0.49, 3.01] | NA | NA | NA | NA | NA | 1.000 | 13003 | 10459 | 0.18 | 0.06 | [0.05, 0.32] | NA | NA | NA | NA | NA | 1.000 | 12931 | 9839 |
| sd(is_B) | 1.37 | 0.23 | [0.99, 1.92] | NA | NA | NA | NA | NA | 1.000 | 25542 | 30110 | 0.38 | 0.06 | [0.28, 0.52] | NA | NA | NA | NA | NA | 1.000 | 16644 | 25124 |
| sd(is_B:pushing_partner_cw) | 1.72 | 0.82 | [0.24, 3.65] | NA | NA | NA | NA | NA | 1.000 | 10498 | 8114 | 0.17 | 0.07 | [0.03, 0.32] | NA | NA | NA | NA | NA | 1.000 | 12628 | 10344 |
| sd(is_B:pushing_self_cw) | 0.85 | 0.33 | [0.29, 1.72] | NA | NA | NA | NA | NA | 1.000 | 14504 | 13114 | 0.05 | 0.04 | [0.00, 0.15] | NA | NA | NA | NA | NA | 1.000 | 19509 | 21091 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.67 | 0.01 | [0.64, 0.70] | NA | NA | NA | NA | NA | 1.000 | 57411 | 29137 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
## [1] 5.75 971.25
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_persuasion", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_obj_log_persuasion, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.84 [1.76, 1.93] 1.36 0.54
## is_B 1.08 [1.05, 1.12] 1.04 0.93
## is_A:persuasion_self_cw 1.13 [1.10, 1.18] 1.06 0.88
## is_B:persuasion_self_cw 1.09 [1.06, 1.13] 1.04 0.92
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 3.53 [3.35, 3.73] 1.88 0.28
## is_A:persuasion_self_cb 2.94 [2.79, 3.10] 1.71 0.34
## is_B:persuasion_self_cb 3.69 [3.50, 3.90] 1.92 0.27
## is_A:persuasion_partner_cb 2.95 [2.80, 3.11] 1.72 0.34
## Tolerance 95% CI
## [0.52, 0.57]
## [0.89, 0.95]
## [0.85, 0.91]
## [0.88, 0.95]
## [0.91, 0.98]
## [0.27, 0.30]
## [0.32, 0.36]
## [0.26, 0.29]
## [0.32, 0.36]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 33, observations = 3343, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001338618 0.005085253
## sample estimates:
## outlier frequency (expected: 0.00303021238408615 )
## 0.009871373
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_persuasion$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_persuasion,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 101.53*** | 6.28 | [ 89.80, 114.88] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 9441 | 16504 |
| is_B | 135.95*** | 7.32 | [122.13, 151.50] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 13358 | 20795 |
| is_A:persuasion_self_cw | 1.05* | 0.02 | [ 1.01, 1.09] | 0.989 | [0.94, 1.07] | 0.839 | 0.083 | Strong Evidence for Null | 1.000 | 36995 | 32045 |
| is_B:persuasion_self_cw | 1.01 | 0.02 | [ 0.97, 1.05] | 0.659 | [0.94, 1.07] | 0.998 | 0.005 | Very Strong Evidence for Null | 1.000 | 35355 | 33024 |
| is_A:persuasion_partner_cw | 1.02 | 0.02 | [ 0.98, 1.06] | 0.792 | [0.94, 1.07] | 0.988 | 0.007 | Very Strong Evidence for Null | 1.000 | 36380 | 30708 |
| is_B:persuasion_partner_cw | 1.03 | 0.02 | [ 0.99, 1.07] | 0.909 | [0.94, 1.07] | 0.969 | 0.014 | Very Strong Evidence for Null | 1.000 | 40799 | 30254 |
| is_A:persuasion_self_cb | 1.25 | 0.21 | [ 0.89, 1.76] | 0.901 | [0.94, 1.07] | 0.128 | 0.033 | Strong Evidence for Null | 1.000 | 10402 | 17765 |
| is_B:persuasion_self_cb | 0.89 | 0.11 | [ 0.70, 1.14] | 0.824 | [0.94, 1.07] | 0.263 | 0.019 | Very Strong Evidence for Null | 1.000 | 14430 | 21279 |
| is_A:persuasion_partner_cb | 0.89 | 0.13 | [ 0.66, 1.20] | 0.787 | [0.94, 1.07] | 0.249 | 0.019 | Very Strong Evidence for Null | 1.000 | 11805 | 19041 |
| is_B:persuasion_partner_cb | 1.29 | 0.19 | [ 0.96, 1.73] | 0.957 | [0.94, 1.07] | 0.081 | 0.062 | Strong Evidence for Null | 1.000 | 13364 | 21848 |
| is_A:day | 0.98 | 0.05 | [ 0.90, 1.08] | 0.628 | [0.94, 1.07] | 0.798 | 0.004 | Very Strong Evidence for Null | 1.000 | 87762 | 28942 |
| is_B:day | 0.95 | 0.04 | [ 0.87, 1.04] | 0.859 | [0.94, 1.07] | 0.612 | 0.007 | Very Strong Evidence for Null | 1.000 | 93845 | 28595 |
| is_A:weartime_self_cw | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 11.432 | Strong Evidence | 1.000 | 79311 | 30606 |
| is_B:weartime_self_cw | 1.00* | 0.00 | [ 1.00, 1.00] | 0.983 | [0.94, 1.07] | 1.000 | 0.039 | Strong Evidence for Null | 1.000 | 80060 | 29963 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.712 | [0.94, 1.07] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.000 | 12377 | 18753 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.584 | [0.94, 1.07] | 1.000 | 0.011 | Very Strong Evidence for Null | 1.000 | 15565 | 23274 |
| sd(is_A) | 0.34 | 0.04 | [0.26, 0.44] | NA | NA | NA | NA | NA | 1.000 | 12417 | 22026 |
| sd(is_B) | 0.28 | 0.04 | [0.22, 0.37] | NA | NA | NA | NA | NA | 1.000 | 14817 | 22063 |
| sd(is_A:persuasion_self_cw) | 0.06 | 0.02 | [0.02, 0.11] | NA | NA | NA | NA | NA | 1.000 | 14879 | 11539 |
| sd(is_B:persuasion_self_cw) | 0.06 | 0.02 | [0.02, 0.11] | NA | NA | NA | NA | NA | 1.000 | 17283 | 14227 |
| sd(is_A:persuasion_partner_cw) | 0.07 | 0.02 | [0.03, 0.12] | NA | NA | NA | NA | NA | 1.000 | 15202 | 15518 |
| sd(is_B:persuasion_partner_cw) | 0.07 | 0.03 | [0.01, 0.13] | NA | NA | NA | NA | NA | 1.000 | 11803 | 13070 |
| sigma | 0.55 | 0.01 | [0.53, 0.56] | NA | NA | NA | NA | NA | 1.000 | 73910 | 28826 |
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pressure", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_obj_log_pressure, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_pressure, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model
## without interaction terms.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.86 [ 1.78, 1.95] 1.36 0.54
## is_B 1.02 [ 1.01, 1.09] 1.01 0.98
## is_A:pressure_self_cw 1.01 [ 1.00, 1.13] 1.01 0.99
## is_B:pressure_self_cw 1.02 [ 1.00, 1.11] 1.01 0.98
## is_A:pressure_partner_cw 1.04 [ 1.01, 1.09] 1.02 0.96
## Tolerance 95% CI
## [0.51, 0.56]
## [0.91, 0.99]
## [0.88, 1.00]
## [0.90, 1.00]
## [0.92, 0.99]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A:pressure_self_cb 8.49 [ 8.01, 9.01] 2.91 0.12
## is_A:pressure_partner_cb 8.81 [ 8.30, 9.34] 2.97 0.11
## Tolerance 95% CI
## [0.11, 0.12]
## [0.11, 0.12]
##
## High Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_B:pressure_partner_cw 10.71 [10.10, 11.37] 3.27 0.09
## is_B:pressure_self_cb 10.70 [10.08, 11.35] 3.27 0.09
## Tolerance 95% CI
## [0.09, 0.10]
## [0.09, 0.10]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:
## 2 of 40000 iterations ended with a divergence (0.005%).
## Try increasing 'adapt_delta' to remove the divergences.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 34, observations = 3337, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001640695 0.004794726
## sample estimates:
## outlier frequency (expected: 0.00308660473479173 )
## 0.01018879
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pressure$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_pressure,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Warning: There were 2 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Model has no intercept. VIFs may not be sensible.
## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model
## without interaction terms.
## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_A 1.86 [ 1.78, 1.95] 1.36 0.54
is_B 1.02 [ 1.01, 1.09] 1.01 0.98
is_A:pressure_self_cw 1.01 [ 1.00, 1.13] 1.01 0.99
is_B:pressure_self_cw 1.02 [ 1.00, 1.11] 1.01 0.98
is_A:pressure_partner_cw 1.04 [ 1.01, 1.09] 1.02 0.96 Tolerance 95% CI [0.51, 0.56] [0.91, 0.99] [0.88, 1.00] [0.90, 1.00] [0.92, 0.99]
Moderate Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_A:pressure_self_cb 8.49 [ 8.01, 9.01] 2.91 0.12
is_A:pressure_partner_cb 8.81 [ 8.30, 9.34] 2.97 0.11 Tolerance 95% CI [0.11, 0.12] [0.11, 0.12]
High Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_B:pressure_partner_cw 10.71 [10.10, 11.37] 3.27 0.09 is_B:pressure_self_cb 10.70 [10.08, 11.35] 3.27 0.09 Tolerance 95% CI [0.09, 0.10] [0.09, 0.10]
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_obj_log_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 101.90*** | 6.39 | [ 89.74, 115.49] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 7210 | 14382 |
| is_B | 137.60*** | 7.65 | [123.21, 153.70] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 10736 | 18368 |
| is_A:pressure_self_cw | 0.99 | 0.05 | [ 0.87, 1.09] | 0.609 | [0.94, 1.07] | 0.764 | 0.006 | Very Strong Evidence for Null | 1.000 | 35375 | 26044 |
| is_B:pressure_self_cw | 0.99 | 0.04 | [ 0.90, 1.08] | 0.616 | [0.94, 1.07] | 0.843 | 0.005 | Very Strong Evidence for Null | 1.000 | 54399 | 31354 |
| is_A:pressure_partner_cw | 1.01 | 0.04 | [ 0.94, 1.09] | 0.603 | [0.94, 1.07] | 0.895 | 0.004 | Very Strong Evidence for Null | 1.000 | 56778 | 29304 |
| is_B:pressure_partner_cw | 1.02 | 0.06 | [ 0.92, 1.17] | 0.667 | [0.94, 1.07] | 0.709 | 0.006 | Very Strong Evidence for Null | 1.000 | 37773 | 26375 |
| is_A:pressure_self_cb | 0.84 | 0.28 | [ 0.43, 1.67] | 0.697 | [0.94, 1.07] | 0.135 | 0.017 | Very Strong Evidence for Null | 1.000 | 17538 | 23698 |
| is_B:pressure_self_cb | 0.94 | 0.22 | [ 0.58, 1.51] | 0.607 | [0.94, 1.07] | 0.210 | 0.012 | Very Strong Evidence for Null | 1.000 | 18333 | 23028 |
| is_A:pressure_partner_cb | 1.25 | 0.33 | [ 0.73, 2.14] | 0.797 | [0.94, 1.07] | 0.136 | 0.022 | Very Strong Evidence for Null | 1.000 | 16840 | 23213 |
| is_B:pressure_partner_cb | 1.20 | 0.36 | [ 0.65, 2.21] | 0.722 | [0.94, 1.07] | 0.139 | 0.014 | Very Strong Evidence for Null | 1.000 | 19355 | 24488 |
| is_A:day | 0.95 | 0.04 | [ 0.87, 1.05] | 0.842 | [0.94, 1.07] | 0.642 | 0.006 | Very Strong Evidence for Null | 1.000 | 83189 | 29167 |
| is_B:day | 0.92 | 0.04 | [ 0.84, 1.01] | 0.961 | [0.94, 1.07] | 0.360 | 0.018 | Very Strong Evidence for Null | 1.000 | 90539 | 28746 |
| is_A:weartime_self_cw | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 13.519 | Strong Evidence | 1.000 | 76844 | 32076 |
| is_B:weartime_self_cw | 1.00 | 0.00 | [ 1.00, 1.00] | 0.973 | [0.94, 1.07] | 1.000 | 0.025 | Very Strong Evidence for Null | 1.000 | 89204 | 27836 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.631 | [0.94, 1.07] | 1.000 | 0.012 | Very Strong Evidence for Null | 1.000 | 9815 | 16497 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.725 | [0.94, 1.07] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.000 | 12846 | 21776 |
| sd(is_A) | 0.34 | 0.04 | [0.27, 0.44] | NA | NA | NA | NA | NA | 1.000 | 10546 | 18623 |
| sd(is_B) | 0.29 | 0.04 | [0.23, 0.38] | NA | NA | NA | NA | NA | 1.000 | 12764 | 19054 |
| sd(is_A:pressure_self_cw) | 0.08 | 0.07 | [0.00, 0.25] | NA | NA | NA | NA | NA | 1.000 | 14433 | 21148 |
| sd(is_B:pressure_self_cw) | 0.05 | 0.04 | [0.00, 0.17] | NA | NA | NA | NA | NA | 1.000 | 21197 | 20973 |
| sd(is_A:pressure_partner_cw) | 0.04 | 0.03 | [0.00, 0.13] | NA | NA | NA | NA | NA | 1.000 | 27025 | 21802 |
| sd(is_B:pressure_partner_cw) | 0.07 | 0.06 | [0.00, 0.27] | NA | NA | NA | NA | NA | 1.000 | 19755 | 22073 |
| sigma | 0.56 | 0.01 | [0.54, 0.57] | NA | NA | NA | NA | NA | 1.000 | 76294 | 29458 |
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_obj_log_pushing, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_pushing, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.00 [1.88, 2.13] 1.41 0.50
## is_B 1.02 [1.00, 1.21] 1.01 0.98
## barriers_self_cw 1.07 [1.03, 1.13] 1.03 0.94
## barriers_self_cb 1.07 [1.04, 1.14] 1.03 0.93
## is_A:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_B:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_A:pushing_partner_cw 1.10 [1.07, 1.17] 1.05 0.91
## is_B:pushing_partner_cw 1.96 [1.85, 2.09] 1.40 0.51
## is_A:pushing_self_cb 1.83 [1.73, 1.94] 1.35 0.55
## is_B:pushing_self_cb 1.96 [1.85, 2.09] 1.40 0.51
## is_A:pushing_partner_cb 1.89 [1.78, 2.01] 1.37 0.53
## Tolerance 95% CI
## [0.47, 0.53]
## [0.83, 1.00]
## [0.88, 0.97]
## [0.88, 0.96]
## [0.89, 0.98]
## [0.89, 0.98]
## [0.86, 0.94]
## [0.48, 0.54]
## [0.51, 0.58]
## [0.48, 0.54]
## [0.50, 0.56]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 13, observations = 1594, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000925345 0.005646173
## sample estimates:
## outlier frequency (expected: 0.00331869510664994 )
## 0.008155583
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pushing$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_pushing,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 106.40*** | 7.11 | [ 93.02, 121.57] | 1.000 | [0.94, 1.06] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 15901 | 24332 |
| is_B | 146.91*** | 9.89 | [128.51, 168.01] | 1.000 | [0.94, 1.06] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 19659 | 25323 |
| barriers_self_cw | 0.95** | 0.02 | [ 0.91, 0.99] | 0.996 | [0.94, 1.06] | 0.647 | 0.177 | Moderate Evidence for Null | 1.000 | 91051 | 28693 |
| barriers_self_cb | 0.82* | 0.08 | [ 0.68, 0.99] | 0.980 | [0.94, 1.06] | 0.070 | 0.115 | Moderate Evidence for Null | 1.000 | 20765 | 25131 |
| is_A:pushing_self_cw | 1.07 | 0.04 | [ 0.99, 1.15] | 0.966 | [0.94, 1.06] | 0.433 | 0.046 | Strong Evidence for Null | 1.000 | 42743 | 30873 |
| is_B:pushing_self_cw | 1.00 | 0.03 | [ 0.94, 1.06] | 0.506 | [0.94, 1.06] | 0.951 | 0.007 | Very Strong Evidence for Null | 1.000 | 42212 | 30965 |
| is_A:pushing_partner_cw | 0.98 | 0.02 | [ 0.94, 1.04] | 0.728 | [0.94, 1.06] | 0.964 | 0.007 | Very Strong Evidence for Null | 1.000 | 67404 | 32270 |
| is_B:pushing_partner_cw | 1.06 | 0.04 | [ 0.99, 1.14] | 0.957 | [0.94, 1.06] | 0.491 | 0.035 | Strong Evidence for Null | 1.000 | 42448 | 31807 |
| is_A:pushing_self_cb | 0.98 | 0.28 | [ 0.55, 1.73] | 0.530 | [0.94, 1.06] | 0.174 | 0.015 | Very Strong Evidence for Null | 1.000 | 16613 | 23372 |
| is_B:pushing_self_cb | 0.98 | 0.31 | [ 0.52, 1.88] | 0.525 | [0.94, 1.06] | 0.154 | 0.018 | Very Strong Evidence for Null | 1.000 | 17222 | 23611 |
| is_A:pushing_partner_cb | 1.16 | 0.36 | [ 0.63, 2.18] | 0.688 | [0.94, 1.06] | 0.140 | 0.018 | Very Strong Evidence for Null | 1.000 | 15659 | 23690 |
| is_B:pushing_partner_cb | 1.42 | 0.41 | [ 0.79, 2.58] | 0.884 | [0.94, 1.06] | 0.080 | 0.033 | Strong Evidence for Null | 1.000 | 17956 | 24013 |
| is_A:day | 1.00 | 0.07 | [ 0.88, 1.14] | 0.523 | [0.94, 1.06] | 0.661 | 0.005 | Very Strong Evidence for Null | 1.000 | 81551 | 28339 |
| is_B:day | 0.89 | 0.06 | [ 0.78, 1.01] | 0.969 | [0.94, 1.06] | 0.186 | 0.031 | Strong Evidence for Null | 1.000 | 80514 | 29652 |
| is_A:weartime_self_cw | 1.00** | 0.00 | [ 1.00, 1.00] | 0.997 | [0.94, 1.06] | 1.000 | 0.237 | Moderate Evidence for Null | 1.000 | 74748 | 33237 |
| is_B:weartime_self_cw | 1.00* | 0.00 | [ 1.00, 1.00] | 0.993 | [0.94, 1.06] | 1.000 | 0.112 | Moderate Evidence for Null | 1.000 | 83919 | 30123 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.703 | [0.94, 1.06] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.000 | 16870 | 23704 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.882 | [0.94, 1.06] | 1.000 | 0.028 | Very Strong Evidence for Null | 1.000 | 19627 | 25505 |
| sd(is_A) | 0.33 | 0.05 | [0.25, 0.44] | NA | NA | NA | NA | NA | 1.000 | 18055 | 25973 |
| sd(is_B) | 0.34 | 0.05 | [0.26, 0.45] | NA | NA | NA | NA | NA | 1.000 | 17389 | 24079 |
| sd(is_A:pushing_self_cw) | 0.10 | 0.05 | [0.01, 0.21] | NA | NA | NA | NA | NA | 1.000 | 9624 | 11024 |
| sd(is_B:pushing_self_cw) | 0.09 | 0.04 | [0.01, 0.17] | NA | NA | NA | NA | NA | 1.000 | 11822 | 9435 |
| sd(is_A:pushing_partner_cw) | 0.03 | 0.03 | [0.00, 0.10] | NA | NA | NA | NA | NA | 1.000 | 23687 | 23076 |
| sd(is_B:pushing_partner_cw) | 0.10 | 0.04 | [0.01, 0.18] | NA | NA | NA | NA | NA | 1.001 | 16684 | 13022 |
| sigma | 0.52 | 0.01 | [0.50, 0.54] | NA | NA | NA | NA | NA | 1.000 | 59441 | 29518 |
## [1] 0 5
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_persuasion", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(mood_gauss_persuasion, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_persuasion, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.17 [1.13, 1.21] 1.08 0.86
## is_B 1.04 [1.02, 1.09] 1.02 0.96
## is_A:persuasion_self_cw 1.06 [1.04, 1.11] 1.03 0.94
## is_B:persuasion_self_cw 1.07 [1.04, 1.11] 1.03 0.93
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 2.85 [2.71, 2.99] 1.69 0.35
## is_A:persuasion_self_cb 2.80 [2.67, 2.94] 1.67 0.36
## Tolerance 95% CI
## [0.82, 0.88]
## [0.92, 0.98]
## [0.90, 0.96]
## [0.90, 0.96]
## [0.92, 0.98]
## [0.33, 0.37]
## [0.34, 0.38]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 47%
## normal 34%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 39, observations = 3688, p-value =
## 2.231e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.007530249 0.014428062
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01057484
summary_mood <- summarize_brms(
mood_gauss_persuasion,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.60*** | 0.12 | [ 3.35, 3.84] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.002 | 3494 | 6838 |
| is_B | 3.88*** | 0.14 | [ 3.60, 4.15] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 4288 | 9553 |
| is_A:persuasion_self_cw | 0.01 | 0.02 | [-0.04, 0.05] | 0.633 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 53463 | 30917 |
| is_B:persuasion_self_cw | -0.01 | 0.02 | [-0.06, 0.03] | 0.717 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 47802 | 29543 |
| is_A:persuasion_partner_cw | 0.05* | 0.02 | [ 0.00, 0.09] | 0.982 | [-0.11, 0.11] | 0.998 | 0.027 | Very Strong Evidence for Null | 1.000 | 45958 | 31085 |
| is_B:persuasion_partner_cw | 0.01 | 0.02 | [-0.04, 0.05] | 0.626 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 55921 | 29328 |
| is_A:persuasion_self_cb | -0.29 | 0.33 | [-0.97, 0.38] | 0.810 | [-0.11, 0.11] | 0.185 | 0.023 | Very Strong Evidence for Null | 1.000 | 5518 | 10646 |
| is_B:persuasion_self_cb | 0.37 | 0.33 | [-0.30, 1.05] | 0.868 | [-0.11, 0.11] | 0.148 | 0.033 | Strong Evidence for Null | 1.001 | 6266 | 11337 |
| is_A:persuasion_partner_cb | 0.41 | 0.29 | [-0.18, 0.99] | 0.916 | [-0.11, 0.11] | 0.120 | 0.041 | Strong Evidence for Null | 1.000 | 5913 | 10851 |
| is_B:persuasion_partner_cb | -0.20 | 0.39 | [-0.98, 0.57] | 0.699 | [-0.11, 0.11] | 0.205 | 0.020 | Very Strong Evidence for Null | 1.000 | 6005 | 11939 |
| is_A:day | 0.45*** | 0.07 | [ 0.30, 0.58] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 51520 | 30176 |
| is_B:day | 0.08 | 0.07 | [-0.06, 0.22] | 0.871 | [-0.11, 0.11] | 0.673 | 0.006 | Very Strong Evidence for Null | 1.000 | 54047 | 28585 |
| sd(is_A) | 0.70 | 0.09 | [0.55, 0.90] | NA | NA | NA | NA | NA | 1.000 | 7045 | 14067 |
| sd(is_B) | 0.79 | 0.10 | [0.63, 1.03] | NA | NA | NA | NA | NA | 1.000 | 7522 | 12132 |
| sd(is_A:pushing_self_cw) | 0.11 | 0.06 | [0.01, 0.23] | NA | NA | NA | NA | NA | 1.000 | 10529 | 10208 |
| sd(is_B:pushing_self_cw) | 0.06 | 0.05 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.000 | 13285 | 16387 |
| sd(is_A:pushing_partner_cw) | 0.13 | 0.05 | [0.04, 0.24] | NA | NA | NA | NA | NA | 1.000 | 17381 | 12874 |
| sd(is_B:pushing_partner_cw) | 0.12 | 0.06 | [0.01, 0.25] | NA | NA | NA | NA | NA | 1.000 | 10329 | 7714 |
| sigma | 0.87 | 0.01 | [0.85, 0.89] | NA | NA | NA | NA | NA | 1.000 | 52078 | 28947 |
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_pressure", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(mood_gauss_pressure, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_pressure, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.19 [1.15, 1.23] 1.09 0.84
## is_B 1.01 [1.00, 1.21] 1.01 0.99
## is_A:pressure_self_cw 1.01 [1.00, 1.12] 1.01 0.99
## is_B:pressure_self_cw 1.01 [1.00, 1.12] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 1.42] 1.00 0.99
## Tolerance 95% CI
## [0.81, 0.87]
## [0.83, 1.00]
## [0.89, 1.00]
## [0.89, 1.00]
## [0.70, 1.00]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_B:pressure_partner_cw 8.40 [7.94, 8.89] 2.90 0.12
## is_A:pressure_self_cb 8.05 [7.61, 8.52] 2.84 0.12
## Tolerance 95% CI
## [0.11, 0.13]
## [0.12, 0.13]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 47%
## normal 34%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 41, observations = 3684, p-value <
## 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.007998021 0.015068026
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01112921
summary_mood <- summarize_brms(
mood_gauss_pressure,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.64*** | 0.13 | [ 3.37, 3.89] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 2966 | 6814 |
| is_B | 3.88*** | 0.14 | [ 3.60, 4.17] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.003 | 3288 | 6776 |
| is_A:pressure_self_cw | -0.06 | 0.06 | [-0.17, 0.06] | 0.835 | [-0.11, 0.11] | 0.845 | 0.005 | Very Strong Evidence for Null | 1.000 | 58461 | 29843 |
| is_B:pressure_self_cw | 0.01 | 0.06 | [-0.10, 0.12] | 0.542 | [-0.11, 0.11] | 0.957 | 0.003 | Very Strong Evidence for Null | 1.000 | 44197 | 29836 |
| is_A:pressure_partner_cw | 0.05 | 0.06 | [-0.06, 0.16] | 0.808 | [-0.11, 0.11] | 0.864 | 0.005 | Very Strong Evidence for Null | 1.000 | 42613 | 30902 |
| is_B:pressure_partner_cw | 0.04 | 0.06 | [-0.07, 0.15] | 0.754 | [-0.11, 0.11] | 0.907 | 0.004 | Very Strong Evidence for Null | 1.000 | 68143 | 31795 |
| is_A:pressure_self_cb | -0.16 | 0.64 | [-1.48, 1.11] | 0.600 | [-0.11, 0.11] | 0.142 | 0.015 | Very Strong Evidence for Null | 1.000 | 10697 | 15995 |
| is_B:pressure_self_cb | 0.12 | 0.57 | [-1.03, 1.31] | 0.588 | [-0.11, 0.11] | 0.157 | 0.017 | Very Strong Evidence for Null | 1.000 | 10971 | 16363 |
| is_A:pressure_partner_cb | 0.22 | 0.50 | [-0.76, 1.27] | 0.674 | [-0.11, 0.11] | 0.164 | 0.016 | Very Strong Evidence for Null | 1.000 | 10831 | 15307 |
| is_B:pressure_partner_cb | -0.12 | 0.72 | [-1.61, 1.35] | 0.566 | [-0.11, 0.11] | 0.126 | 0.016 | Very Strong Evidence for Null | 1.000 | 11345 | 16597 |
| is_A:day | 0.41*** | 0.07 | [ 0.28, 0.55] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 49901 | 30936 |
| is_B:day | 0.09 | 0.07 | [-0.05, 0.22] | 0.897 | [-0.11, 0.11] | 0.647 | 0.006 | Very Strong Evidence for Null | 1.000 | 50514 | 29576 |
| sd(is_A) | 0.71 | 0.09 | [0.57, 0.93] | NA | NA | NA | NA | NA | 1.001 | 5363 | 10000 |
| sd(is_B) | 0.81 | 0.10 | [0.65, 1.06] | NA | NA | NA | NA | NA | 1.001 | 6354 | 13168 |
| sd(is_A:pushing_self_cw) | 0.12 | 0.06 | [0.01, 0.25] | NA | NA | NA | NA | NA | 1.001 | 9570 | 7098 |
| sd(is_B:pushing_self_cw) | 0.07 | 0.05 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.000 | 12893 | 14758 |
| sd(is_A:pushing_partner_cw) | 0.13 | 0.05 | [0.04, 0.24] | NA | NA | NA | NA | NA | 1.000 | 13798 | 9709 |
| sd(is_B:pushing_partner_cw) | 0.12 | 0.06 | [0.02, 0.25] | NA | NA | NA | NA | NA | 1.000 | 12399 | 9266 |
| sigma | 0.87 | 0.01 | [0.85, 0.89] | NA | NA | NA | NA | NA | 1.000 | 43311 | 29442 |
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_pushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(mood_gauss_pushing, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_pushing, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.35 [1.29, 1.43] 1.16 0.74
## is_B 1.01 [1.00, 1.48] 1.01 0.99
## barriers_self_cw 1.09 [1.05, 1.15] 1.04 0.92
## barriers_self_cb 1.09 [1.05, 1.15] 1.04 0.92
## is_A:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.95
## is_B:pushing_self_cw 1.08 [1.05, 1.14] 1.04 0.93
## is_A:pushing_partner_cw 1.09 [1.05, 1.15] 1.04 0.92
## is_B:pushing_partner_cw 1.50 [1.42, 1.58] 1.22 0.67
## is_A:pushing_self_cb 1.54 [1.46, 1.62] 1.24 0.65
## Tolerance 95% CI
## [0.70, 0.77]
## [0.68, 1.00]
## [0.87, 0.95]
## [0.87, 0.95]
## [0.89, 0.98]
## [0.88, 0.96]
## [0.87, 0.95]
## [0.63, 0.70]
## [0.62, 0.68]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 41%
## normal 41%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 11, observations = 1776, p-value =
## 0.00112
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.003095804 0.011055146
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.006193694
summary_mood <- summarize_brms(
mood_gauss_pushing,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.66*** | 0.12 | [ 3.41, 3.90] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 7565 | 14933 |
| is_B | 3.84*** | 0.14 | [ 3.55, 4.13] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 7731 | 13693 |
| barriers_self_cw | -0.29*** | 0.03 | [-0.35, -0.23] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 78719 | 29655 |
| barriers_self_cb | -0.50* | 0.22 | [-0.95, -0.05] | 0.984 | [-0.11, 0.11] | 0.039 | 0.167 | Moderate Evidence for Null | 1.000 | 8107 | 14716 |
| is_A:pushing_self_cw | 0.09 | 0.05 | [ 0.00, 0.18] | 0.969 | [-0.11, 0.11] | 0.659 | 0.034 | Strong Evidence for Null | 1.000 | 42568 | 31360 |
| is_B:pushing_self_cw | -0.01 | 0.04 | [-0.09, 0.06] | 0.640 | [-0.11, 0.11] | 0.990 | 0.005 | Very Strong Evidence for Null | 1.000 | 48493 | 32295 |
| is_A:pushing_partner_cw | 0.09 | 0.04 | [ 0.00, 0.18] | 0.974 | [-0.11, 0.11] | 0.678 | 0.039 | Strong Evidence for Null | 1.000 | 32188 | 29785 |
| is_B:pushing_partner_cw | 0.13* | 0.05 | [ 0.02, 0.23] | 0.985 | [-0.11, 0.11] | 0.376 | 0.077 | Strong Evidence for Null | 1.000 | 33217 | 29086 |
| is_A:pushing_self_cb | -0.23 | 0.48 | [-1.20, 0.73] | 0.681 | [-0.11, 0.11] | 0.162 | 0.015 | Very Strong Evidence for Null | 1.000 | 11236 | 17666 |
| is_B:pushing_self_cb | 0.90 | 0.67 | [-0.47, 2.27] | 0.905 | [-0.11, 0.11] | 0.053 | 0.047 | Strong Evidence for Null | 1.000 | 9401 | 14271 |
| is_A:pushing_partner_cb | 0.77 | 0.52 | [-0.29, 1.83] | 0.927 | [-0.11, 0.11] | 0.056 | 0.042 | Strong Evidence for Null | 1.001 | 9141 | 15831 |
| is_B:pushing_partner_cb | -0.32 | 0.63 | [-1.57, 0.96] | 0.693 | [-0.11, 0.11] | 0.124 | 0.022 | Very Strong Evidence for Null | 1.000 | 11283 | 17409 |
| is_A:day | 0.34*** | 0.10 | [ 0.15, 0.54] | 1.000 | [-0.11, 0.11] | 0.008 | 1.941 | Weak Evidence | 1.000 | 70757 | 30332 |
| is_B:day | 0.08 | 0.10 | [-0.11, 0.27] | 0.793 | [-0.11, 0.11] | 0.609 | 0.006 | Very Strong Evidence for Null | 1.000 | 69006 | 28830 |
| sd(is_A) | 0.64 | 0.09 | [0.50, 0.84] | NA | NA | NA | NA | NA | 1.001 | 10062 | 18475 |
| sd(is_B) | 0.79 | 0.10 | [0.63, 1.04] | NA | NA | NA | NA | NA | 1.000 | 10080 | 17533 |
| sd(is_A:pushing_self_cw) | 0.09 | 0.07 | [0.01, 0.24] | NA | NA | NA | NA | NA | 1.000 | 12455 | 15086 |
| sd(is_B:pushing_self_cw) | 0.07 | 0.05 | [0.00, 0.20] | NA | NA | NA | NA | NA | 1.000 | 14139 | 16967 |
| sd(is_A:pushing_partner_cw) | 0.12 | 0.06 | [0.01, 0.24] | NA | NA | NA | NA | NA | 1.000 | 13704 | 10047 |
| sd(is_B:pushing_partner_cw) | 0.14 | 0.07 | [0.02, 0.29] | NA | NA | NA | NA | NA | 1.000 | 14270 | 12569 |
| sigma | 0.83 | 0.01 | [0.80, 0.86] | NA | NA | NA | NA | NA | 1.000 | 60197 | 29169 |
## [1] 0 5
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_persuasion", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(is_reactance_persuasion)
DHARMa.check_brms.all(is_reactance_persuasion, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.82 [1.71, 1.96] 1.35 0.55
## is_B 1.39 [1.32, 1.49] 1.18 0.72
## is_A:persuasion_self_cw 1.15 [1.10, 1.23] 1.07 0.87
## is_B:persuasion_self_cw 1.04 [1.01, 1.14] 1.02 0.96
## is_A:persuasion_partner_cw 1.06 [1.02, 1.14] 1.03 0.95
## is_B:persuasion_partner_cw 1.51 [1.43, 1.62] 1.23 0.66
## is_A:persuasion_self_cb 1.78 [1.67, 1.91] 1.33 0.56
## Tolerance 95% CI
## [0.51, 0.58]
## [0.67, 0.76]
## [0.82, 0.91]
## [0.88, 0.99]
## [0.88, 0.98]
## [0.62, 0.70]
## [0.52, 0.60]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 44%
## cauchy 12%
## gamma 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 81%
## beta-binomial 12%
## binomial 6%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.001322751
summary_is_reactance <- summarize_brms(
is_reactance_persuasion,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.18*** | 0.08 | [0.07, 0.43] | 1.000 | [0.83, 1.20] | 0.000 | 36.028 | Very Strong Evidence | 1.000 | 24897 | 29836 |
| is_B | 0.30** | 0.12 | [0.12, 0.65] | 0.999 | [0.83, 1.20] | 0.004 | 2.802 | Weak Evidence | 1.000 | 29961 | 29498 |
| is_A:persuasion_self_cw | 0.90 | 0.13 | [0.64, 1.18] | 0.773 | [0.83, 1.20] | 0.669 | 0.071 | Strong Evidence for Null | 1.000 | 22239 | 26651 |
| is_B:persuasion_self_cw | 0.78 | 0.12 | [0.54, 1.02] | 0.965 | [0.83, 1.20] | 0.310 | 0.293 | Moderate Evidence for Null | 1.000 | 23382 | 19058 |
| is_A:persuasion_partner_cw | 0.89 | 0.16 | [0.58, 1.29] | 0.740 | [0.83, 1.20] | 0.579 | 0.077 | Strong Evidence for Null | 1.000 | 30521 | 22671 |
| is_B:persuasion_partner_cw | 1.27 | 0.24 | [0.86, 1.97] | 0.897 | [0.83, 1.20] | 0.356 | 0.165 | Moderate Evidence for Null | 1.000 | 25245 | 25126 |
| is_A:persuasion_self_cb | 5.13 | 4.34 | [0.96, 35.39] | 0.972 | [0.83, 1.20] | 0.025 | 0.569 | Weak Evidence for Null | 1.000 | 18611 | 24445 |
| is_B:persuasion_self_cb | 2.76 | 2.14 | [0.62, 15.25] | 0.910 | [0.83, 1.20] | 0.078 | 0.216 | Moderate Evidence for Null | 1.000 | 20515 | 24894 |
| is_A:persuasion_partner_cb | 1.58 | 1.13 | [0.40, 7.92] | 0.746 | [0.83, 1.20] | 0.169 | 0.114 | Moderate Evidence for Null | 1.000 | 16564 | 19332 |
| is_B:persuasion_partner_cb | 2.13 | 1.94 | [0.34, 14.06] | 0.798 | [0.83, 1.20] | 0.109 | 0.122 | Moderate Evidence for Null | 1.000 | 18086 | 24106 |
| is_A:day | 2.57 | 1.36 | [0.91, 7.43] | 0.962 | [0.83, 1.20] | 0.058 | 0.194 | Moderate Evidence for Null | 1.000 | 56235 | 31190 |
| is_B:day | 1.04 | 0.59 | [0.33, 3.17] | 0.523 | [0.83, 1.20] | 0.248 | 0.045 | Strong Evidence for Null | 1.000 | 52331 | 30861 |
| sd(is_A) | 1.38 | 0.42 | [0.68, 2.40] | NA | NA | NA | NA | NA | 1.000 | 14059 | 19287 |
| sd(is_B) | 1.36 | 0.33 | [0.82, 2.18] | NA | NA | NA | NA | NA | 1.000 | 16662 | 23957 |
| sd(is_A:persuasion_self_cw) | 0.31 | 0.18 | [0.03, 0.73] | NA | NA | NA | NA | NA | 1.001 | 9510 | 11805 |
| sd(is_B:persuasion_self_cw) | 0.35 | 0.26 | [0.02, 0.95] | NA | NA | NA | NA | NA | 1.000 | 5934 | 13700 |
| sd(is_A:persuasion_partner_cw) | 0.47 | 0.27 | [0.05, 1.19] | NA | NA | NA | NA | NA | 1.000 | 12291 | 11597 |
| sd(is_B:persuasion_partner_cw) | 0.65 | 0.29 | [0.19, 1.40] | NA | NA | NA | NA | NA | 1.000 | 16040 | 15521 |
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_pressure", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(is_reactance_pressure)
DHARMa.check_brms.all(is_reactance_pressure, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.83 [1.72, 1.96] 1.35 0.55
## is_B 1.02 [1.00, 1.23] 1.01 0.98
## is_A:pressure_self_cw 1.01 [1.00, 1.59] 1.01 0.99
## is_B:pressure_self_cw 1.01 [1.00, 2.01] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 36.14] 1.00 0.99
## is_B:pressure_partner_cw 2.51 [2.33, 2.71] 1.58 0.40
## is_A:pressure_self_cb 2.97 [2.76, 3.22] 1.72 0.34
## Tolerance 95% CI
## [0.51, 0.58]
## [0.82, 1.00]
## [0.63, 1.00]
## [0.50, 1.00]
## [0.03, 1.00]
## [0.37, 0.43]
## [0.31, 0.36]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 34%
## beta 16%
## cauchy 12%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 81%
## beta-binomial 12%
## binomial 6%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 23 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.001322751
summary_is_reactance <- summarize_brms(
is_reactance_pressure,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.15*** | 0.06 | [0.06, 0.32] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 27201 | 29563 |
| is_B | 0.21*** | 0.07 | [0.10, 0.40] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 33935 | 30484 |
| is_A:pressure_self_cw | 2.48 | 1.27 | [0.82, 11.13] | 0.957 | [0.83, 1.20] | 0.046 | 0.724 | Weak Evidence for Null | 1.000 | 21167 | 16301 |
| is_B:pressure_self_cw | 1.84 | 0.67 | [0.78, 4.37] | 0.935 | [0.83, 1.20] | 0.096 | 0.359 | Weak Evidence for Null | 1.000 | 27629 | 24426 |
| is_A:pressure_partner_cw | 0.87 | 0.65 | [0.09, 3.62] | 0.581 | [0.83, 1.20] | 0.205 | 0.128 | Moderate Evidence for Null | 1.000 | 24646 | 22567 |
| is_B:pressure_partner_cw | 1.56 | 0.77 | [0.39, 5.39] | 0.805 | [0.83, 1.20] | 0.170 | 0.128 | Moderate Evidence for Null | 1.000 | 30905 | 28417 |
| is_A:pressure_self_cb | 16.91 | 28.48 | [0.78, 1025.11] | 0.965 | [0.83, 1.20] | 0.018 | 0.537 | Weak Evidence for Null | 1.000 | 26287 | 20773 |
| is_B:pressure_self_cb | 23.56* | 34.63 | [1.93, 874.04] | 0.994 | [0.83, 1.20] | 0.006 | 1.817 | Weak Evidence | 1.000 | 25130 | 23908 |
| is_A:pressure_partner_cb | 1.05 | 1.44 | [0.05, 17.51] | 0.513 | [0.83, 1.20] | 0.105 | 0.120 | Moderate Evidence for Null | 1.000 | 24236 | 22816 |
| is_B:pressure_partner_cb | 0.36 | 0.65 | [0.00, 9.70] | 0.725 | [0.83, 1.20] | 0.071 | 0.105 | Moderate Evidence for Null | 1.000 | 29760 | 25453 |
| is_A:day | 2.24 | 1.18 | [0.79, 6.37] | 0.936 | [0.83, 1.20] | 0.087 | 0.138 | Moderate Evidence for Null | 1.000 | 59228 | 29462 |
| is_B:day | 1.34 | 0.70 | [0.48, 3.72] | 0.710 | [0.83, 1.20] | 0.234 | 0.050 | Strong Evidence for Null | 1.000 | 64117 | 29915 |
| sd(is_A) | 1.46 | 0.39 | [0.84, 2.44] | NA | NA | NA | NA | NA | 1.000 | 14625 | 22279 |
| sd(is_B) | 1.14 | 0.28 | [0.67, 1.83] | NA | NA | NA | NA | NA | 1.000 | 15965 | 24630 |
| sd(is_A:pressure_self_cw) | 1.30 | 1.07 | [0.06, 3.91] | NA | NA | NA | NA | NA | 1.000 | 9827 | 16955 |
| sd(is_B:pressure_self_cw) | 1.03 | 0.55 | [0.17, 2.56] | NA | NA | NA | NA | NA | 1.000 | 13715 | 12493 |
| sd(is_A:pressure_partner_cw) | 1.81 | 1.16 | [0.13, 4.53] | NA | NA | NA | NA | NA | 1.000 | 13418 | 11212 |
| sd(is_B:pressure_partner_cw) | 0.74 | 0.66 | [0.03, 2.89] | NA | NA | NA | NA | NA | 1.000 | 22546 | 21566 |
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_pushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(is_reactance_pushing)
DHARMa.check_brms.all(is_reactance_pushing, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.85 [1.75, 1.97] 1.36 0.54
## is_B 1.06 [1.03, 1.13] 1.03 0.94
## barriers_self_cw 1.08 [1.05, 1.15] 1.04 0.92
## barriers_self_cb 1.24 [1.19, 1.31] 1.11 0.81
## is_A:pushing_self_cw 1.14 [1.10, 1.21] 1.07 0.87
## is_B:pushing_self_cw 1.01 [1.00, 1.89] 1.00 0.99
## is_A:pushing_partner_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_B:pushing_partner_cw 1.53 [1.45, 1.62] 1.24 0.65
## is_A:pushing_self_cb 1.33 [1.27, 1.40] 1.15 0.75
## Tolerance 95% CI
## [0.51, 0.57]
## [0.89, 0.97]
## [0.87, 0.96]
## [0.77, 0.84]
## [0.83, 0.91]
## [0.53, 1.00]
## [0.89, 0.98]
## [0.62, 0.69]
## [0.71, 0.79]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 44%
## cauchy 12%
## tweedie 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 75%
## beta-binomial 16%
## binomial 9%
##
## Divergences:
## 0 of 40000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
##
## DHARMa outlier test based on exact binomial test with
## approximate expectations
##
## data: model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.000000000 0.007561553
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0
summary_is_reactance <- summarize_brms(
is_reactance_pushing,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.13*** | 0.07 | [0.04, 0.33] | 1.000 | [0.83, 1.20] | 0.000 | 75.098 | Very Strong Evidence | 1.000 | 21167 | 22208 |
| is_B | 0.16*** | 0.07 | [0.06, 0.35] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 28688 | 28514 |
| barriers_self_cw | 0.80 | 0.17 | [0.53, 1.19] | 0.865 | [0.83, 1.20] | 0.388 | 0.098 | Strong Evidence for Null | 1.000 | 43333 | 29118 |
| barriers_self_cb | 2.74 | 2.06 | [0.61, 12.70] | 0.910 | [0.83, 1.20] | 0.077 | 0.221 | Moderate Evidence for Null | 1.000 | 20027 | 24196 |
| is_A:pushing_self_cw | 1.35 | 0.25 | [0.93, 2.00] | 0.944 | [0.83, 1.20] | 0.258 | 0.217 | Moderate Evidence for Null | 1.000 | 35407 | 26611 |
| is_B:pushing_self_cw | 1.44 | 0.29 | [0.96, 2.23] | 0.964 | [0.83, 1.20] | 0.167 | 0.440 | Weak Evidence for Null | 1.000 | 25829 | 25680 |
| is_A:pushing_partner_cw | 0.77 | 0.24 | [0.32, 1.36] | 0.815 | [0.83, 1.20] | 0.338 | 0.129 | Moderate Evidence for Null | 1.000 | 19974 | 15679 |
| is_B:pushing_partner_cw | 1.15 | 0.23 | [0.75, 1.74] | 0.752 | [0.83, 1.20] | 0.517 | 0.079 | Strong Evidence for Null | 1.000 | 35235 | 26529 |
| is_A:pushing_self_cb | 38.45* | 65.62 | [2.09, 2544.92] | 0.993 | [0.83, 1.20] | 0.005 | 2.260 | Weak Evidence | 1.000 | 12493 | 16241 |
| is_B:pushing_self_cb | 2.37 | 3.65 | [0.13, 75.77] | 0.717 | [0.83, 1.20] | 0.082 | 0.123 | Moderate Evidence for Null | 1.000 | 17952 | 21560 |
| is_A:pushing_partner_cb | 0.08 | 0.16 | [0.00, 4.46] | 0.898 | [0.83, 1.20] | 0.032 | 0.268 | Moderate Evidence for Null | 1.000 | 16536 | 21879 |
| is_B:pushing_partner_cb | 1.25 | 1.87 | [0.05, 26.76] | 0.561 | [0.83, 1.20] | 0.097 | 0.100 | Strong Evidence for Null | 1.000 | 19752 | 22818 |
| is_A:day | 1.28 | 0.94 | [0.29, 5.36] | 0.633 | [0.83, 1.20] | 0.184 | 0.059 | Strong Evidence for Null | 1.000 | 39476 | 29465 |
| is_B:day | 1.40 | 0.92 | [0.39, 5.18] | 0.696 | [0.83, 1.20] | 0.194 | 0.064 | Strong Evidence for Null | 1.000 | 47667 | 31034 |
| sd(is_A) | 1.52 | 0.52 | [0.61, 2.78] | NA | NA | NA | NA | NA | 1.000 | 8591 | 8207 |
| sd(is_B) | 1.26 | 0.38 | [0.61, 2.16] | NA | NA | NA | NA | NA | 1.000 | 14530 | 19907 |
| sd(is_A:pushing_self_cw) | 0.26 | 0.23 | [0.01, 0.91] | NA | NA | NA | NA | NA | 1.000 | 10496 | 17794 |
| sd(is_B:pushing_self_cw) | 0.56 | 0.30 | [0.07, 1.32] | NA | NA | NA | NA | NA | 1.001 | 8769 | 8991 |
| sd(is_A:pushing_partner_cw) | 0.64 | 0.46 | [0.04, 1.94] | NA | NA | NA | NA | NA | 1.000 | 10726 | 13385 |
| sd(is_B:pushing_partner_cw) | 0.28 | 0.25 | [0.01, 1.05] | NA | NA | NA | NA | NA | 1.000 | 14535 | 18263 |
Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)